home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 247_01 / hilbert.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-04-17  |  2.4 KB  |  110 lines
  1. /*
  2.  * Solve set of linear equations involving
  3.  * a Hilbert matrix
  4.  * i.e. solves   Hx=b, where b is the vector [1,1,1....1]
  5.  */
  6.  
  7. #include <stdio.h>
  8. #include "miracl.h"
  9.  
  10. flash A[50][50];
  11. flash b[50];
  12.  
  13. bool gauss(A,b,n)
  14. flash A[][50];
  15. flash b[];
  16. int n;
  17. { /* solve Ax=b using Gaussian elimination *
  18.    * solution x returned in b              */
  19.     int i,j,k,m;
  20.     bool ok;
  21.     flash w,s;
  22.     w=mirvar(0);
  23.     s=mirvar(0);
  24.     ok=TRUE;
  25.     for (i=0;i<n;i++)
  26.         copy(b[i],A[i][n]);
  27.     for (i=0;i<n;i++)
  28.     { /* Gaussian elimination */
  29.         m=i;
  30.         for (j=i+1;j<n;j++)
  31.         {
  32.             absol(A[j][i],w);
  33.             absol(A[m][i],s);
  34.             if (fcomp(w,s)>0) m=j;
  35.         }
  36.         if (m!=i) for (k=i;k<=n;k++)
  37.         {
  38.             copy(A[i][k],w);
  39.             copy(A[m][k],A[i][k]);
  40.             copy(w,A[m][k]);
  41.         }
  42.         if (size(A[i][i])==0)
  43.         {
  44.             ok=FALSE;
  45.             break;
  46.         }
  47.         for (j=i+1;j<n;j++)
  48.         {
  49.             fdiv(A[j][i],A[i][i],s);
  50.             for (k=n;k>=i;k--)
  51.             {
  52.                 fmul(s,A[i][k],w);
  53.                 fsub(A[j][k],w,A[j][k]);
  54.             }
  55.         }
  56.     }
  57.     if (ok) for (j=n-1;j>=0;j--)
  58.     { /* Backward substitution */
  59.         zero(s);
  60.         for (k=j+1;k<n;k++)
  61.         {
  62.             fmul(A[j][k],b[k],w);
  63.             fadd(s,w,s);
  64.         }
  65.         fsub(A[j][n],s,w);
  66.         if (size(A[j][j])==0)
  67.         {
  68.             ok=FALSE;
  69.             break;
  70.         } 
  71.         fdiv(w,A[j][j],b[j]);
  72.     }
  73.     free(s);
  74.     free(w);
  75.     return ok;
  76. }
  77.  
  78. main()
  79. { /* solve set of linear equations */
  80.     int i,j,n;
  81.     mirsys(20,MAXBASE);
  82.     do
  83.     {
  84.         printf("Order of Hilbert matrix H= ");
  85.         scanf("%d",&n);
  86.         getchar();
  87.     } while (n<2 || n>49);
  88.     for (i=0;i<n;i++)
  89.     {
  90.         A[i][n]=mirvar(0);
  91.         b[i]=mirvar(1);
  92.         for (j=0;j<n;j++)
  93.         {
  94.             A[i][j]=mirvar(0);
  95.             fconv(1,i+j+1,A[i][j]);
  96.         }
  97.     }
  98.     if (gauss(A,b,n))
  99.     {
  100.         printf("\nSolution is\n");
  101.         for (i=0;i<n;i++)
  102.         {
  103.             printf("x[%d] = ",i+1);
  104.             cotnum(b[i],stdout);
  105.         }
  106.         if (EXACT) printf("Result is exact!\n");
  107.     }
  108.     else printf("H is singular!\n");
  109. }
  110.